/***********************************************
CHILD OUTCOMES

THIS DO-FILE ANALYSES ALL CHILD LEVEL FINAL OUTCOMES
- ASSESSMENTS OF CHILD DEVELOPMENT
- ANTHROPOMETRICS
***********************************************/


clear all
set more off
capture log close

 if c(username) == "alison_a" {
        global ROOT "C:\Users\alison_a\Dropbox\HI_shared (1)\"
		global ROOT1 "C:\Users\alison_a\Dropbox\HI_shared\"
    }

if c(username) == "sonya_k" {
        global ROOT "C:\Users\sonya_k\Dropbox\HI_shared\"
		global ROOT1 "C:\Users\sonya_k\Dropbox\HI_private\"
    }
	
if c(username) == "Sonya" {
        global ROOT "C:\Users\Sonya\Dropbox\HI_shared"
		global ROOT1 "C:\Users\Sonya\Dropbox\HI_private\"
    }

if c(username) == "sofyakrutikova" {
        global ROOT "/Users/sofyakrutikova/Dropbox/HI_shared/endline report/"
		global ROOT1 "/Users/sofyakrutikova/Dropbox/HI_private/"
    }

			
global output 				"$ROOT/endline report/draft 2"
global do 					"$ROOT/do-files"
global data_constructed 	"C:\Users\alison_a\Dropbox\HI_shared (1)\endline report\FINAL\DATA"
global rw 					"$ROOT/data/RW/T-value distributions"
		
		
do "$do\stata programs for drawing tables\programs for endline tables, standardisations etc"		

	* load data
	use "$data_constructed\child_analysis_data"
	
	* rename
	rename mac_say_and_und_0* macarthur_0*

	
	* set up common heterogeneity variables
	*************************************************
	*************************************************
	* which variables
	*1 - age
	sum agedays_0, detail
	gen older=(agedays_0>=r(p50)) if agedays_0~=.
	tab older	

	*2-gender - already defined
	
	*3-developmental level at BL
	* develop factor for development at BL and split in half
	factor asq_com_0_as_z asq_fm_0_as_z asq_mg_0_as_z asq_ps_0_as_z asq_soc_0_as_z macarthur_0_as_z // one factor
	predict temp
	sum temp, detail
	gen bl_dev=(temp>=r(p50)) if temp~=.
	tab bl_dev
	drop temp
	
	*4-maternal education at BL
	sum educmadre_0, detail
	gen educ_moth=(educmadre_0>=12) if educmadre~=.
	tab educ_moth
	
	*5 - Level of Stimulation in the Home at Baseline
	factor fci_*_0 // use main factor
	predict temp
	sum temp, detail
	gen bl_stim=(temp>=r(p50)) if temp~=.
	tab bl_stim
	drop temp
	
	* 6- size of childcare centre
	sum num_children_0, detail
	gen bl_size=(num_children_0>=r(p50)) if num_children_0~=.
	tab bl_size
	
	* 7- ICBF standards - BL
	sum icbf_qual_0, detail
	gen bl_qual=(icbf_qual_0>=r(p50)) if icbf_qual_0~=.
	tab bl_qual
	
	* 8 - teach education 
	sum educmaestra_0, detail
	gen bl_teach_educ=(educmaestra_0>=r(p50)) if educmaestra_0~=.
	tab bl_teach_educ
	
	*9- teacher experience
	sum expmaestra_0, detail
	gen bl_teach_exp=(expmaestra_0>=r(p50)) if expmaestra_0~=.
	tab bl_teach_exp
	
	*10- teacher burnout
	sum total_burnout_0, detail
	gen bl_teach_burn=(total_burnout_0>=r(p50)) if total_burnout_0~=.
	tab bl_teach_burn
	
	*11 - teacher job satisfaction
	sum jobsatisfaction_0, detail
	gen bl_teach_satis=(jobsatisfaction_0>=r(p50)) if jobsatisfaction_0~=.
	tab bl_teach_satis
	
	*12 - teacher depression	
	sum depresion_0, detail
	gen bl_teach_dep=(depresion_0>=r(p50)) if depresion_0~=.
	tab bl_teach_dep
	
	*13 - productive routines	
	sum  rutinas_P1, detail
	gen bl_prod_rout=( rutinas_P1>=r(p50)) if  rutinas_P1~=.
	tab bl_prod_rout
	
	
	*define common control variables for all child level regressions
	**************************************************************************
	tab dane, gen(dane_)
	
	* control variable globals
	global bl_tests "asq_com_0_as_z asq_mg_0_as_z asq_ps_0_as_z asq_soc_0_as_z asq_fm_0_as_z macarthur_0_as_z" // could change to as !!!!!!!!!!!!!!!!!!!!!!!
	global bl_tests_temp "asq_com_0_ats_z asq_mg_0_ats_z asq_ps_0_ats_z asq_soc_0_ats_z asq_fm_0_ats_z macarthur_0_ats_z" // could change to as !!!!!!!!!!!!!!!!!!!!!!!
	global controls_child "age_months_ex*t male dane_*"
	global controls_child_tests "male dane_*" 
	
		
	* interviewer dummies
	replace cc_ent=80727225 if cc_ent==0
	tab cc_ent, gen(cc_ent_)

	tab cc_tester, gen(cc_tester_)
	
		
	* CHILD DEVELOPMENT OUTCOMES
	*****************************
	*****************************
	
	*NUTRITIONAL STATUS
	****************************************************************************
	
	* z-scores
	************
	global vars "zwei zlen zbmi zwfl"

	* define labels
	local zwei 		"Weight for age z-score"
	local zlen 		"Length/Height for age z-score"
	local zbmi		"BMI for age z-score"
	local zwfl		"Weight for length z-score"

	* label standardised vars
	foreach x in $vars {
		label var `x' "``x''"
		}
	
	* create 2 sided test table --- no correction for multiple testing
	*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	myttests $vars, itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_ent_*   $controls_child) bl(y)
	table_2s, filename(zscores)
	
	* create sidak corrected pvalues
	*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	tempfile sidak_temp
	tempname sidak
	postfile `sidak' str22 test itt2_pval_2s  itt3_pval_2s   itt23_pval_2s  using `sidak_temp', replace

	foreach x in zwei zlen zbmi zwfl { 
		reg `x' itt2 itt3 `x'_0 $controls_child  dane_* cc_ent_* , cluster(ll_inst )
		test itt2=0
		local itt2_p2s=r(p)
		
		test itt3=0
		local itt3_p2s=r(p)
				
		test itt2=itt3	
		local itt23_p2s=r(p)
		
		post `sidak' ("`x'")  (`itt2_p2s')  (`itt3_p2s')  (`itt23_p2s') 
		}
		
		postclose `sidak'

	tempfile temp3
	save `temp3'
	
	use `sidak_temp', clear
	foreach i of numlist 2 3 23 {
			qqvalue itt`i'_pval_2s, method(sidak) qvalue(itt`i'_sidak_2s)
			}
		tempfile sidak
		save `sidak', replace
			
	use `temp3', clear
	myttests  zwei zlen zbmi zwfl, itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_ent_*   $controls_child) bl(y)
	preserve
	use `sidak', clear
	mySIDAKpvalues_return zwei zlen zbmi zwfl
	restore
	table_RW_2s, filename(SIDAK_anthro)
	

	* HETEROGENEITY of impacts
	******************
		
	global hvars "older male bl_dev educ_moth bl_*"
	global vars "zwei zlen zbmi zwfl"
	
	 *setup
	hetero_setup $vars, hetero(older) label0("Younger children") label1("Older children")
	hetero_setup $vars, hetero(male) label0("Female") label1("Male")
	hetero_setup $vars, hetero(bl_dev) label0("Less developed") label1("More developed")
	hetero_setup $vars, hetero(educ_moth) label0("Mother <12y educ") label1("Mother >=12y educ")
	hetero_setup $vars, hetero(bl_stim) label0("Low home stimulation") label1("High home stimulation")
	hetero_setup $vars, hetero(bl_size) label0("Smaller centres") label1("Larger centres")
	hetero_setup $vars, hetero(bl_qual) label0("Lower ICBF quality score") label1("Higher ICBF quality score")
	hetero_setup $vars, hetero(bl_teach_educ) label0("Teacher less educ") label1("Teacher more educ")
	hetero_setup $vars, hetero(bl_teach_exp) label0("Teacher less exp") label1("Teacher more exp")
	hetero_setup $vars, hetero(bl_teach_burn) label0("Lower teacher burnout score") label1("Higher teacher burnout score")
	hetero_setup $vars, hetero(bl_teach_satis) label0("Lower teacher job satis score") label1("Higher teacher job satis score")
	hetero_setup $vars, hetero(bl_teach_dep) label0("Lower teacher depression score") label1("Higher teacher depression score")
	hetero_setup $vars, hetero(bl_prod_rout) label0("Fewer productive routines") label1("More productive routines")

	* define labels
	local zwei 		"Weight for age z-score"
	local zlen 		"Length/Height for age z-score"
	local zbmi		"BMI for age z-score"
	local zwfl		"Weight for length z-score"
	
	* relabel to make bold font
	foreach x in $vars {
		label var `x' "\textbf{``x''}"
		}
		
	* analysis
	hetero $vars , bl(y) hvars($hvars) itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_ent_* $controls_child) sided(2)
	table_het_2s, filename(zscores_het)


			
	* binary nutritional indicators
	***********************************
	
	* drop nutritional indicators with under 10% or over 90% prevalence
	drop nutr_vlow_wfh nutr_acut* nutr_obes* nutr_glob* // only one case
	
	* unadjusted for multiple hypothesis testing
	myttests_bin nutr_chron_mal-nutr_ov_bmi , itt2(itt2) itt3(itt3) bl(y) controls(cc_ent_* $controls_child) cluster(ll_inst  )
	table_bin_2s, filename(nutr_measures)
	
	* create Sidak p-values for multiple testing
	cap drop nutr_vlow_wfh nutr_acut* nutr_obes* nutr_glob* // only one case

	tempfile sidak_temp
	tempname sidak
	postfile `sidak' str22 test itt2_pval_2s  itt3_pval_2s   itt23_pval_2s  using `sidak_temp', replace
	
	foreach x of varlist nutr_chron_mal-nutr_ov_bmi { 
		logit `x'  itt2 itt3 `x'_0 cc_ent_* $controls_child  dane_*, cluster(ll_inst) 
		test itt2=0
		local itt2_p2s=r(p)
		test itt3=0
		local itt3_p2s=r(p)		
		test itt2=itt3	
		local itt23_p2s=r(p)
		post `sidak' ("`x'")  (`itt2_p2s')  (`itt3_p2s')  (`itt23_p2s') 
		}

	postclose `sidak'

	tempfile temp3
	save `temp3'
	
	use `sidak_temp', clear
	foreach i of numlist 2 3 23 {
			qqvalue itt`i'_pval_2s, method(sidak) qvalue(itt`i'_sidak_2s)
			}
	tempfile sidak
	save `sidak', replace
			
	use `temp3', clear
	myttests_bin  nutr_chron_mal nutr_low_height nutr_ad_height nutr_ad_weight ///
		nutr_risklow_wfh  nutr_ad_wfh   nutr_overweight nutr_ov_bmi, itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_ent_*   $controls_child) bl(y)
	
	preserve
	use `sidak', clear
	mySIDAKpvalues_return nutr_chron_mal nutr_low_height nutr_ad_height nutr_ad_weight ///
		nutr_risklow_wfh  nutr_ad_wfh   nutr_overweight nutr_ov_bmi
	restore
	
	table_bin_RW_2s, filename(SIDAK_nutr_measures)
	
	
	* ALL CHILD ASSESSMENTS
	****************************************************************************
	****************************************************************************
	
	* CHILD COGNITIVE AND LITERACY DEVELOPMENT
	******************************************
	******************************************
	tempfile keep
	save `keep' // before dropping outliers
	
	* drop outliers - zscore of <-3 in any one of cognitive tests -- not on 5 and 21 - we will drop these
	
	foreach x in IRT_tvip_as_z  wm12_as_z  IRT_wm14_as_z  IRT_wm17_as_z   ptt_as_z  IRT_dab_as_z  { // CHANGED TO IRT -- 11.12.15
		di "`x'"
		drop if `x'<-3
		count if `x'>3 & `x'~=.
		}

	
	* ALL TESTS
	********************
	* LABEL TESTS
	*********************
	*label var wm5_as_z "Woodcock Munoz - 5"
	label var wm12_as_z "Woodcock Munoz - 12"
	label var IRT_wm14_as_z "Woodcock Munoz - 14"
	label var IRT_wm17_as_z "Woodcock Munoz - 17"
	label var IRT_wm21_as_z "Woodcock Munoz - 21"
	label var IRT_tvip_as_z "TVIP"
	label var ptt_as_z "Pencil Tapping Test"
	label var IRT_dab_as_z "Daberon - total"	
	
	
	* TREATMENT EFFECTS
	**********************
	* NB. OTHER SPECIFICATIONS (INCLUDING THAT USED IN FIRST REPORT) IN APPENDIX OF THE DO FILE
	* HERE JUST INCLUDE FULL IRT SPECIFICATION (FOR ALL ASSESSMENTS WHERE IRT IS APPLICABLE)
	* TESTER EFFECTS REMOVED VIA TESTER DUMMIES ON THE 
		
	****************************************************************************
	* full IRT 
	****************************************************************************
	
	* COG FACTOR 1
	factor IRT_tvip_as_z IRT_dab_as_z ptt_as_z wm12_as_z IRT_wm14_as_z IRT_wm17_as_z, mineigen(1) blanks(0.25)
	predict factor_cog1_as_fullirt
	lab var factor_cog1_as_fullirt "Cog - inc PTT - full irt - alt controls"	
	sum factor_cog1_as_fullirt if itt2==0 & itt3==0
	replace factor_cog1_as_fullirt=(factor_cog1_as_fullirt-r(mean))/r(sd)
	
	* export factor loadings
	mat E=e(Ev)
	estadd scalar ev=E[1,1]
	esttab using "$output/cog_fact_loadings1.tex",  scalars("ev Eigenvalue") sfmt(3) replace gaps nolines cells("L[1](transpose fmt(3)) ") collabel("Factor Loading" )   noobs nonumber nomtitle label
	
	label var factor_cog1 "Cog, Lang and Sch (all measures, limited sample)"
	kdensity factor_cog1, xtitle("Factor value") ytitle("Density") graphregion(color(white)) title("") // title("Cognitive Development, Language Development and School Readiness (all measures) (factor)")
	graph export "$output/graph_cog_factor1.emf", replace
	
		
	* COG FACTOR 2 - EXC PTT
	factor IRT_tvip_as_z IRT_dab_as_z  wm12_as_z IRT_wm14_as_z IRT_wm17_as_z, mineigen(1) blanks(0.25)
	predict factor_cog2_as_fullirt
	label var factor_cog2 "Cog, Lang and Sch (exc. PTT, full sample)"
	sum factor_cog2_as_fullirt if itt2==0 & itt3==0
	replace factor_cog2_as_fullirt=(factor_cog2_as_fullirt-r(mean))/r(sd)	
	
	* export factor loadings
	mat E=e(Ev)
	estadd scalar ev=E[1,1]
	esttab using "$output/cog_fact_loadings2.tex",  scalars("ev Eigenvalue") sfmt(3) replace gaps nolines cells("L[1](transpose fmt(3)) ") collabel("Factor Loading"  )   noobs nonumber nomtitle label
	
	kdensity factor_cog2, xtitle("Standardised Score") ytitle("Density") graphregion(color(white)) title("")
	graph export "$output/graph_cog_factor2.emf", replace
	
	
	* COG FACTOR 3 - LIMITED SAMPLE
	gen factor_cog3_as_fullirt=factor_cog2_as_fullirt if factor_cog1_as_fullirt~=.
	lab var factor_cog3_as_fullirt "Cog, Lang and Sch (exc. PTT, limited sample)"	
	
	sum factor_cog3_as_fullirt if itt2==0 & itt3==0
	replace factor_cog3_as_fullirt=(factor_cog3_as_fullirt-r(mean))/r(sd)
	
	
	* PRELIT
	factor IRT_tvip_as_z   IRT_wm14_as_z IRT_wm17_as_z, mineigen(1) blanks(0.25)
	predict factor_prelit_as_fullirt
	label var factor_prelit "Pre-literacy skills"

	sum factor_prelit_as_fullirt if itt2==0 & itt3==0
	replace factor_prelit_as_fullirt=(factor_prelit_as_fullirt-r(mean))/r(sd)	

	* export factor loadings
	mat E=e(Ev)
	estadd scalar ev=E[1,1]
	esttab using "$output/prelit_fact_loadings.tex",  scalars("ev Eigenvalue") sfmt(3) replace gaps nolines cells("L[1](transpose fmt(3)) ") collabel("Factor Loading" )   noobs nonumber nomtitle label
	
	kdensity factor_prelit, xtitle("Pre-literacy skills") ytitle("Density") graphregion(color(white)) title("") 
	graph export "$output/prelit_cog_factor.emf", replace
	
	
	* DRAW TABLE
	****************
	****************
	myttests factor_cog1_as_fullirt  factor_cog3_as_fullirt factor_cog2_as_fullirt factor_prelit_as_fullirt  ,  ///
		itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_tester_* $bl_tests  $controls_child_tests)
	table_std_1s, filename(FACTORS_IRT_TESTERDUMMIES)
	
	
	
	***************************************************************************
	* INDIVIDUAL TESTS ( irt scores where applicable and tester effects as dummies) - no adjustmet for multiple testing
	******************************
	******************************
	myttests IRT_tvip_as_z  IRT_wm14_as_z IRT_wm17_as_z wm12_as_z ptt_as_z IRT_dab_as_z  ,  ///
		itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_tester_* $bl_tests  $controls_child_tests)
	table_std_1s, filename(noRW_TESTS_IRT_TESTERDUMMIES)
	
	
	* RW multiple hypothesis testing
	******************************************
	tempfile keep2
	save `keep2'
		
	/*
	myRWpvalues_setup_alt IRT_tvip_as_z  IRT_wm14_as_z IRT_wm17_as_z wm12_as_z ptt_as_z IRT_dab_as_z , ///
			b(5000) sided(1) itt2(itt2) itt3(itt3) cluster(ll_inst )  ///
			controls(cc_tester_* $bl_tests  $controls_child_tests) setname(child_tests_1201)
	*/		
	
	use `keep2'		

	* RW CORRECTED pvalues TABLE FOR IMPACTS ON CHILD DEVELOPMENT ASSESSMENTS
	myttests  IRT_tvip_as_z  wm12_as_z IRT_wm14_as_z IRT_wm17_as_z  ptt_as_z IRT_dab_as_z, itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_tester_* $bl_tests  $controls_child_tests)
	myRWpvalues_return IRT_tvip_as_z  IRT_wm14_as_z IRT_wm17_as_z wm12_as_z ptt_as_z IRT_dab_as_z, setname(child_tests_1201)
	table_std_RW_1s, filename(RW_TESTS_IRT_TESTERDUMMIES)
	
	
	
	* DABERON SUBSCALES
	****************************************************************************
	*label daberon sub-scales
	label var dab_as_z   				"Daberon total"
	label var dab_bodyparts_as_z		"Daberon - bodyparts"
	label var dab_colour_as_z			"Daberon - colour concepts"
	label var dab_number_as_z			"Daberon - number concepts"	
    label var dab_prep_as_z				"Daberon - prepositions"
	label var dab_directions_as_z		"Daberon - following directions"
	label var dab_knowledge_as_z		"Daberon - general knowledge"
	label var dab_plurals_as_z			"Daberon - plurals"
	label var dab_visual_as_z			"Daberon - visual perception"
	label var dab_categories_as_z		"Daberon - visual perception"
 
	
	myttests  dab*as_z, itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_tester_* $bl_tests  $controls_child_tests)
	table_std_1s, filename(daberon)
	
	* RW corrected pvalues
	/* 	
	myRWpvalues_setup_alt dab*as_z , ///
			b(5000) sided(1) itt2(itt2) itt3(itt3) cluster(ll_inst )  ///
			controls(cc_tester_* $bl_tests  $controls_child_tests) setname(daberon_1201)
			
	
	myttests  dab*as_z, itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_tester_* $bl_tests  $controls_child_tests)
	myRWpvalues_return dab*as_z, setname(daberon_1201)
	table_std_1s, filename(RW_daberon)
	
	*/
	
	use `keep2'	
	
	
	* COEFFICIENT PLOT BY AGE
	****************************************************************************
			
	* in 100 day age bands
	* split into THREE bands and 
	
	preserve
	centile agemons_0 , centile(0 33) 
	keep if agemons_0<=r(c_2) & agemons_0>r(c_1)
	local minage=r(c_1)
	local maxage=r(c_2)
	regress factor_cog2 itt* cc_tester_* $bl_tests  $controls_child_tests , cluster(ll_inst) level(90)	
	regsave itt3, ci double addlabel(min_age,`minage', max_age, `maxage', factor, 2) level(90)	
	tempfile q1
	save `q1'
	restore
	

	preserve
	centile agemons_0 , centile(33 66) 
	keep if agemons_0<=r(c_2) & agemons_0>r(c_1)
	local minage=r(c_1)
	local maxage=r(c_2)
	regress factor_cog2 itt* cc_tester_* $bl_tests  $controls_child_tests , cluster(ll_inst) level(90)	
	regsave itt3, ci double addlabel(min_age,`minage', max_age, `maxage', factor, 2) level(90)	
	tempfile q2
	save `q2'
	restore
	
	preserve
	centile agemons_0 , centile(66 100) 
	keep if agemons_0<=r(c_2) & agemons_0>r(c_1)
	local minage=r(c_1)
	local maxage=r(c_2)
	regress factor_cog2 itt* cc_tester_* $bl_tests  $controls_child_tests, cluster(ll_inst) level(90)	
	regsave itt3, ci double addlabel(min_age,`minage', max_age, `maxage', factor, 2) level(90)	
	tempfile q3
	save `q3'
	restore

	

	preserve
	centile agemons_0 , centile(33 66) 
	keep if agemons_0<=r(c_2) & agemons_0>r(c_1)
	local minage=r(c_1)
	local maxage=r(c_2)
	regress factor_cog1 itt* cc_tester_* $bl_tests  $controls_child_tests , cluster(ll_inst) level(90)	
	regsave itt3, ci double addlabel(min_age,`minage', max_age, `maxage', factor, 1) level(90)	
	tempfile q2_1
	save `q2_1'
	restore
	
	preserve
	centile agemons_0 , centile(66 100) 
	keep if agemons_0<=r(c_2) & agemons_0>r(c_1)
	local minage=r(c_1)
	local maxage=r(c_2)
	regress factor_cog1 itt* cc_tester_* $bl_tests  $controls_child_tests, cluster(ll_inst) level(90)	
	regsave itt3, ci double addlabel(min_age,`minage', max_age, `maxage', factor, 1) level(90)	
	tempfile q3_1
	save `q3_1'
	restore
	
	/*preserve
	centile agedays_0 , centile(60 80) 
	keep if agedays_0<=r(c_2) & agedays_0>r(c_1)
	local minage=r(c_1)
	local maxage=r(c_2)
	regress factor_cog2 itt* $controls_child $bl_tests cc_ent_*, cluster(ll_inst) level(90)	
	regsave itt3, ci double addlabel(min_age,`minage', max_age, `maxage')
	tempfile q4
	save `q4'
	restore
	
	preserve
	centile agedays_0 , centile(80 100) 
	keep if agedays_0<=r(c_2) & agedays_0>r(c_1)
	local minage=r(c_1)
	local maxage=r(c_2)
	regress factor_cog2 itt* $controls_child $bl_tests cc_ent_*, cluster(ll_inst) level(90)	
	regsave itt3, ci double addlabel(min_age,`minage', max_age, `maxage')
	tempfile q5
	save `q5'
	restore	*/
	
	preserve
	use `q1', clear
	append using `q2'
	append using `q3'
	*append using `q1_2'
	append using `q2_1'
	append using `q3_1'
	
	gen mid_age=(min_age+max_age)/2
	twoway (scatter coef mid_age if factor==1, xscale(r(20,40)) yscale(r(-0.125 0.3)) xlabel(20(5)40)) (rcap ci_lower	ci_upper  mid_age if factor==1, xscale(r(20 40))), xtitle("Age in months at baseline") ytitle("Estimated impact of HIM+FE") legend(lab(1 "Point estimate") lab(2 "90% Confidence Interval")) graphregion(color(white)) 
	graph export "$output/ate_coeffplot_1.emf", replace

	twoway (scatter coef mid_age if factor==2, xscale(r(20,40)) yscale(r(-0.125 0.3)) xlabel(20(5)40))  (rcap ci_lower	ci_upper  mid_age if factor==2), xtitle("Age in months at baseline") ytitle("Estimated impact of HIM+FE") legend(lab(1 "Point estimate") lab(2 "90% Confidence Interval")) graphregion(color(white)) 
	graph export "$output/ate_coeffplot_2.emf", replace
	
	restore
	

	
	* HETEROGENEITY
	******************
	
	* heterogeneity by compliance 
	******************************
	merge m:1 ll_inst using "$data_constructed\Compliance data"
	
	* HIM compliance
	************************
	* define more compliant and less compliant HIM centers - those that met at least half the required ratio for each of the three categories
	*
	gen more_compliant_him=(staffing_se_ratio_fte>=0.0025 & staffing_ne_ratio_fte>=0.0025 & staffing_ta_ratio_fte>=0.01 & (itt2==1 | itt3==1 ))
	tab more_compliant_him
	
	gen itt2_morecomp_him=more_compliant_him*itt2
	gen itt2_lesscomp_him=(more_compliant_him==0)*itt2
	
	gen itt3_morecomp_him=more_compliant_him*itt3
	gen itt3_lesscomp_him=(more_compliant_him==0)*itt3

	reg factor_cog1  itt*comp*him cc_tester_* $controls_child $bl_tests, cluster(ll_inst)
	esttab using "$output/het by HIM compliance - fac1.csv", replace nogaps
	* one sided p-values
	foreach coef in itt2_morecomp itt2_lesscomp itt3_morecomp  itt3_lesscomp {
		qui test `coef'
		local sign = sign(_b[`coef'_him])
		display "Ho: `coef' <= 0  p-value = " ttail(r(df_r),`sign'*sqrt(r(F)))
		*display "Ho: coef >= 0  p-value = " 1-ttail(r(df_r),`sign'*sqrt(r(F)))
		}
	test itt2_morecomp=itt2_lesscomp	
	test itt3_morecomp=itt3_lesscomp
	
	reg factor_cog2  itt*comp*him cc_tester_* $controls_child $bl_tests, cluster(ll_inst)
	esttab using "$output/het by HIM compliance - fac2.csv", replace nogaps
	* one sided p-values
	foreach coef in itt2_morecomp itt2_lesscomp itt3_morecomp  itt3_lesscomp {
		qui test `coef'
		local sign = sign(_b[`coef'_him])
		display "Ho: `coef' <= 0  p-value = " ttail(r(df_r),`sign'*sqrt(r(F)))
		*display "Ho: coef >= 0  p-value = " 1-ttail(r(df_r),`sign'*sqrt(r(F)))
		}
	test itt2_morecomp=itt2_lesscomp	
	test itt3_morecomp=itt3_lesscomp	

	reg factor_cog3  itt*comp*him cc_tester_* $controls_child $bl_tests, cluster(ll_inst)
	esttab using "$output/het by HIM compliance - fac3.csv", replace nogaps
	* one sided p-values
	foreach coef in itt2_morecomp itt2_lesscomp itt3_morecomp  itt3_lesscomp {
		qui test `coef'
		local sign = sign(_b[`coef'_him])
		display "Ho: `coef' <= 0  p-value = " ttail(r(df_r),`sign'*sqrt(r(F)))
		*display "Ho: coef >= 0  p-value = " 1-ttail(r(df_r),`sign'*sqrt(r(F)))
		}
	test itt2_morecomp=itt2_lesscomp	
	test itt3_morecomp=itt3_lesscomp	
	
	reg factor_prelit  itt*comp*him cc_tester_* $controls_child $bl_tests, cluster(ll_inst)
	esttab using "$output/het by HIM compliance - fac prelit.csv", replace nogaps
	* one sided p-values
	foreach coef in itt2_morecomp itt2_lesscomp itt3_morecomp  itt3_lesscomp {
		qui test `coef'
		local sign = sign(_b[`coef'_him])
		display "Ho: `coef' <= 0  p-value = " ttail(r(df_r),`sign'*sqrt(r(F)))
		*display "Ho: coef >= 0  p-value = " 1-ttail(r(df_r),`sign'*sqrt(r(F)))
		}
	test itt2_morecomp=itt2_lesscomp	
	test itt3_morecomp=itt3_lesscomp	


	* FE compliance
	**********************
	* define more compliant and less compliant FE centers - if more than a third of teachers received FE training
	***************	
	gen more_compliant_fe=( ratio_fe_training_teachers>=0.33 & itt3==1)
	tab more_compliant_fe if itt3==1

	
	gen itt3_morecomp_fe=more_compliant_fe*itt3
	gen itt3_lesscomp_fe=(more_compliant_fe==0)*itt3	

	reg factor_cog1  itt*comp*fe itt2 cc_tester_* $controls_child $bl_tests, cluster(ll_inst)
	esttab using "$output/het by fe compliance - fac1.csv", replace nogaps
	* one sided p-values
	foreach coef in   itt3_morecomp  itt3_lesscomp {
		qui test `coef'
		local sign = sign(_b[`coef'_fe])
		display "Ho: `coef' <= 0  p-value = " ttail(r(df_r),`sign'*sqrt(r(F)))
		*display "Ho: coef >= 0  p-value = " 1-ttail(r(df_r),`sign'*sqrt(r(F)))
		}
	test itt3_morecomp=itt3_lesscomp
	
	reg factor_cog2  itt*comp*fe itt2 cc_tester_* $controls_child $bl_tests, cluster(ll_inst)
	esttab using "$output/het by fe compliance - fac2.csv", replace nogaps
	* one sided p-values
	foreach coef in   itt3_morecomp  itt3_lesscomp {
		qui test `coef'
		local sign = sign(_b[`coef'_fe])
		display "Ho: `coef' <= 0  p-value = " ttail(r(df_r),`sign'*sqrt(r(F)))
		*display "Ho: coef >= 0  p-value = " 1-ttail(r(df_r),`sign'*sqrt(r(F)))
		}
	test itt3_morecomp=itt3_lesscomp	

	reg factor_cog3  itt*comp*fe itt2 cc_tester_* $controls_child $bl_tests, cluster(ll_inst)
	esttab using "$output/het by fe compliance - fac3.csv", replace nogaps
	* one sided p-values
	foreach coef in   itt3_morecomp  itt3_lesscomp {
		qui test `coef'
		local sign = sign(_b[`coef'_fe])
		display "Ho: `coef' <= 0  p-value = " ttail(r(df_r),`sign'*sqrt(r(F)))
		*display "Ho: coef >= 0  p-value = " 1-ttail(r(df_r),`sign'*sqrt(r(F)))
		}
	test itt3_morecomp=itt3_lesscomp	
	
	reg factor_prelit  itt*comp*fe itt2 cc_tester_* $controls_child $bl_tests, cluster(ll_inst)
	esttab using "$output/het by fe compliance - fac prelit.csv", replace nogaps
	* one sided p-values
	foreach coef in   itt3_morecomp  itt3_lesscomp {
		qui test `coef'
		local sign = sign(_b[`coef'_fe])
		display "Ho: `coef' <= 0  p-value = " ttail(r(df_r),`sign'*sqrt(r(F)))
		*display "Ho: coef >= 0  p-value = " 1-ttail(r(df_r),`sign'*sqrt(r(F)))
		}
	test itt3_morecomp=itt3_lesscomp
	

	* heterogeneity by child, teacher and center characteristics
	*************************************************************
	rename factor_prelit_as_fullirt factor_prelit  
	rename factor_cog1_as_fullirt factor_cog1  
	rename factor_cog2_as_fullirt factor_cog2  	
	
	global hvars "older male bl_dev educ_moth bl_stim " // 
	global vars "factor_cog1 factor_cog2 factor_prelit IRT_tvip_as_z  IRT_wm14_as_z IRT_wm17_as_z wm12_as_z ptt_as_z IRT_dab_as_z" //
		
	*setup
	hetero_setup $vars, hetero(older) label0("Younger children") label1("Older children")
	hetero_setup $vars, hetero(male) label0("Female") label1("Male")
	hetero_setup $vars, hetero(bl_dev) label0("Less developed") label1("More developed")
	hetero_setup $vars, hetero(educ_moth) label0("Mother <12y educ") label1("Mother >=12y educ")
	hetero_setup $vars, hetero(bl_stim) label0("Low home stimulation") label1("High home stimulation")
	hetero_setup $vars, hetero(bl_size) label0("Smaller centres") label1("Larger centres")
	hetero_setup $vars, hetero(bl_qual) label0("Lower ICBF quality score") label1("Higher ICBF quality score")
	hetero_setup $vars, hetero(bl_teach_educ) label0("Teacher less educ") label1("Teacher more educ")
	hetero_setup $vars, hetero(bl_teach_exp) label0("Teacher less exp") label1("Teacher more exp")
	hetero_setup $vars, hetero(bl_teach_burn) label0("Lower teacher burnout score") label1("Higher teacher burnout score")
	hetero_setup $vars, hetero(bl_teach_satis) label0("Lower teacher job satis score") label1("Higher teacher job satis score")
	hetero_setup $vars, hetero(bl_teach_dep) label0("Lower teacher depression score") label1("Higher teacher depression score")
	hetero_setup $vars, hetero(bl_prod_rout) label0("Fewer productive routines") label1("More productive routines")

	* analysis
	set trace off
	hetero factor_cog1 ,  hvars(older male  educ_moth bl_*) itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_tester_* $controls_child $bl_tests) sided(1)
	table_het_1s, filename(cog_het_1)
	
	hetero factor_cog2 ,  hvars(older male  educ_moth bl_*) itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_tester_* $controls_child $bl_tests) sided(1)
	table_het_1s, filename(cog_het_2)

			
	/* heterogeneity  tabs for other factors - not included in report
	hetero factor_cog2 ,  hvars( older male bl_*) itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_tester_* $controls_child $bl_tests) sided(1)
	table_het_1s, filename(cog_het_2)

	hetero factor_prelit ,  hvars( older male bl_*) itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_tester_* $controls_child $bl_tests) sided(1)
	table_het_1s, filename(prelit_het)
	*/
	
	* heterogeneity on all tests by age
	hetero IRT_tvip_as_z  IRT_wm14_as_z IRT_wm17_as_z wm12_as_z ptt_as_z IRT_dab_as_z,  hvars( older) itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_tester_* $controls_child $bl_tests) sided(1)
	table_het_1s, filename(all_tests_het_age)
	
	rename factor_prelit factor_prelit_as_fullirt 
	rename factor_cog1 factor_cog1_as_fullirt 
	rename factor_cog2 factor_cog2_as_fullirt 
	
	
	
	* ASQ-SE
	***********
	rename asqse*0_as_z asqse*as_z_0
	rename asqse_adapt_as_z_0  asqse_adapt_funct_as_z_0 
	rename asqse_int_as_z_0 asqse_interaction_as_z_0   
	
	label var asqse_self_reg_as_z "Self Regulation"
	label var asqse_compliance_as_z "Compliance"
	label var asqse_commun_as_z  "Communication"
	label var asqse_adapt_funct_as_z  "Adaptive Functioning"
	label var asqse_autonomy_as_z "Autonomy"
	label var asqse_affect_as_z   "Affect"
	label var asqse_interaction_as_z  "Interaction with People"
	label var asqse_as_z  "\textbf{Total}"
	
		* socio-emotional 	
	
	estpost correlate asqse*_as_z , matrix listwise
	esttab using "$output/correlation_matrix_tests_se.tex", replace nolines unstack not noobs  label  nonumbers nomtitles
	esttab using "$output/correlation_matrix_tests_se.rtf", replace nolines unstack not noobs  label  nonumbers nomtitles
	
	
	* table - uncorrected for multiple testing
	myttests asqse*_as_z, itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_ent_* $controls_child) bl(y)
	table_std_1s, filename(asq_se)
	
		
	* MULTIPLE HYPOTHESIS TESTING - SIDAK 
	*******************************
	tempfile sidak_temp
	tempname sidak
		postfile `sidak' str22 test itt2_pval_2s  itt3_pval_2s   itt23_pval_2s  using `sidak_temp', replace

	*one sided pvalues	
	foreach x of varlist asqse*_as_z { 
		reg `x' itt2 itt3 `x'_0 $controls_child  dane_* cc_ent_* , cluster(ll_inst )
		
		test itt2=0
		local sign_itt2=sign(_b[itt2])
				if r(F)~=. {
					local itt2_p2s	=ttail(r(df_r),`sign_itt2'*sqrt(r(F)))
					}
				if r(F)==. {
					local itt2_p2s	= 1-normal(`sign_itt2'*sqrt(r(chi2)))
					}				

		test itt3=0
		local sign_itt3=sign(_b[itt3])
				if r(F)~=. {
					local itt3_p2s	=ttail(r(df_r),`sign_itt2'*sqrt(r(F)))
					}
				if r(F)==. {
					local itt3_p2s	= 1-normal(`sign_itt2'*sqrt(r(chi2)))
					}				

		test itt2=itt3
		local sign_itt23=sign(_b[itt3]-_b[itt2])
				if r(F)~=. {
					local itt23_p2s	=ttail(r(df_r),`sign_itt2'*sqrt(r(F)))
					}
				if r(F)==. {
					local itt23_p2s	= 1-normal(`sign_itt2'*sqrt(r(chi2)))
					}
					
		
		post `sidak' ("`x'")  (`itt2_p2s')  (`itt3_p2s')  (`itt23_p2s') 
		}
		
		postclose `sidak'

	tempfile temp3
	save `temp3'
	
	use `sidak_temp', clear
	foreach i of numlist 2 3 23 {
			*qqvalue itt`i'_pval_2s, method(holm) qvalue(itt`i'_holm_2s)
			qqvalue itt`i'_pval_2s, method(sidak) qvalue(itt`i'_sidak_2s)
			}
		tempfile sidak
		
		save `sidak', replace
			
	use `temp3', clear
	myttests  asqse*_as_z, itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_ent_*   $controls_child) bl(y)
	preserve
	use `sidak', clear
	mySIDAKpvalues_return asqse_self_reg_as_z asqse_compliance_as_z asqse_commun_as_z  ///
		asqse_adapt_funct_as_z asqse_autonomy_as_z  asqse_affect_as_z  asqse_interaction_as_z  ///
		asqse_as_z  
	restore
	table_std_RW_1s, filename(SIDAK_asqse)
	
		
	* heterogeneity
	****************
	global vars "asqse_as_z " //
	label var asqse_as_z "ASQ Socio Emotional"
	
	*setup
	hetero_setup $vars, hetero(older) label0("Younger children") label1("Older children")
	hetero_setup $vars, hetero(male) label0("Female") label1("Male")
	hetero_setup $vars, hetero(bl_dev) label0("Less developed") label1("More developed")
	hetero_setup $vars, hetero(educ_moth) label0("Mother <12y educ") label1("Mother >=12y educ")
	hetero_setup $vars, hetero(bl_stim) label0("Low home stimulation") label1("High home stimulation")
	hetero_setup $vars, hetero(bl_size) label0("Smaller centres") label1("Larger centres")
	hetero_setup $vars, hetero(bl_qual) label0("Lower ICBF quality score") label1("Higher ICBF quality score")
	hetero_setup $vars, hetero(bl_teach_educ) label0("Teacher less educ") label1("Teacher more educ")
	hetero_setup $vars, hetero(bl_teach_exp) label0("Teacher less exp") label1("Teacher more exp")
	hetero_setup $vars, hetero(bl_teach_burn) label0("Lower teacher burnout score") label1("Higher teacher burnout score")
	hetero_setup $vars, hetero(bl_teach_satis) label0("Lower teacher job satis score") label1("Higher teacher job satis score")
	hetero_setup $vars, hetero(bl_teach_dep) label0("Lower teacher depression score") label1("Higher teacher depression score")
	hetero_setup $vars, hetero(bl_prod_rout) label0("Fewer productive routines") label1("More productive routines")

	* analysis
	set trace off
	hetero asqse_as_z ,  hvars(older male  educ_moth bl_*) itt2(itt2) itt3(itt3) cluster(ll_inst ) controls(cc_ent_* $controls_child $bl_tests) sided(1)
	table_het_1s, filename(asqse_het_1)
	

